The first time you use a new version of R you will have to install all your packages. Try installing the raster and rgdal packages if you don’t have them already:

if(!require("raster")) install.packages("raster")
## Loading required package: raster
## Loading required package: sp
if(!require("rgdal")) install.packages("rgdal")
## Loading required package: rgdal
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-3, (SVN revision 1187)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /home/mmann1123/.local/share/proj:/usr/share/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
if(!require("lwgeom")) install.packages("lwgeom")
## Loading required package: lwgeom
## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.10.2, PROJ 8.2.1
if(!require("sf")) install.packages("sf")
## Loading required package: sf
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
if(!require("dplyr")) install.packages("dplyr")
## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
if(!require("osmdata")) install.packages("osmdata")
## Loading required package: osmdata
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

Raster data is stored in a format similar to that of a matrix of numbers, except you have added data about the coordinate references system (projection).

Let’s read in one raster and take a look

Examine Raster Data

Notice that the working directory doesn’t need to be a path directly below your files of interest:

CWD_apr = raster('./data/raster/CWD/cwd2010apr.tif')

Get info about the raster object

CWD_apr
## class      : RasterLayer 
## dimensions : 1120, 872, 976640  (nrow, ncol, ncell)
## resolution : 1080, 1080  (x, y)
## extent     : -375034, 566726, -616284.9, 593315.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : cwd2010apr.tif 
## names      : cwd2010apr 
## values     : 0, 150.9694  (min, max)

Plot the image:

plot(CWD_apr)

there are a variety of different plot options. Here are a few important ones

plot(CWD_apr,main = 'Title')

plot(CWD_apr,
     col= heat.colors(10), 
     main='Water Deficit' ) # 10 designates the number of hex color codes that get generated

Now use ?plot to change col to a different preset color scheme

?plot
## Help on topic 'plot' was found in the following packages:
## 
##   Package               Library
##   terra                 /home/mmann1123/R/x86_64-pc-linux-gnu-library/4.1
##   raster                /home/mmann1123/R/x86_64-pc-linux-gnu-library/4.1
##   graphics              /usr/lib/R/library
##   base                  /usr/lib/R/library
##   sf                    /home/mmann1123/R/x86_64-pc-linux-gnu-library/4.1
## 
## 
## Using the first match ...

Raster Math

Rasters like matrices can make basic math functions very easy. Note that rasters can be used in basic operations like CWD_apr *2 or CWD_apr^2.

However teo rasters with different extents, cell sizes, or projections CAN’T be used. You would need to reprojected and aligned first, this can be done using the projectRaster() function, run ?projectRaster for help.

CWD_apr_sq = CWD_apr^2
plot(CWD_apr_sq, main='CWD Squared')

Now lets add two different rasters. Use the raster() function to read in the CWD for May

Now let’s add two months of data values:

#CWD_may = raster('fill/in/the/path/file.tif')

CWD_A_M = CWD_apr + CWD_may
plot(CWD_A_M)

Let’s reproject that raster into a UTM zone 10 and write it back out

# Define the Proj.4 spatial reference 
# https://spatialreference.org/ref/epsg/4326/
sr = "EPSG:4326" 

# Project Raster
projected_raster = raster::projectRaster(CWD_apr, crs = sr)

# Write the RasterLayer to disk (See datatype documentation for other formats)
writeRaster(projected_raster, filename="./output/CWD_apr_UTM.tif",  overwrite=TRUE)

Raster Values

We can also directly set the cell values of rasters. For instance, it is often useful to have a raster filled with a fixed value, such as 0. In this example we will use the setValues() function to update the values of CWD_may to all be 0, and store it in a variable called zeros.

zeros = setValues(CWD_may,0)
plot(zeros)

Note that zeros raster has the same extent, resolution, and projection as CWD_may

print(zeros)
## class      : RasterLayer 
## dimensions : 1120, 872, 976640  (nrow, ncol, ncell)
## resolution : 1080, 1080  (x, y)
## extent     : -375034, 566726, -616284.9, 593315.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : cwd2010may 
## values     : 0, 0  (min, max)
print(CWD_may)
## class      : RasterLayer 
## dimensions : 1120, 872, 976640  (nrow, ncol, ncell)
## resolution : 1080, 1080  (x, y)
## extent     : -375034, 566726, -616284.9, 593315.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : cwd2010may.tif 
## names      : cwd2010may

We can also use boolean queries to update certain values. Here we update certain values of CWD_may to 5000 where CWD_may>150 is TRUE with the following syntax:

CWD_may[CWD_may>150] = 5000
plot(CWD_may)

Note that we have now changed the values of CWD_may… if we want to original values again, we will need to read them back in again with CWD_may = raster('./data/raster/CWD/cwd2010may.tif').

Challenge Question

In this section we will try to get the average water deficit for California for the whole year of 2010. To do this we will need to isolate the files we are interested in using - those ending with .tif. We will use the function list.files() which creates a list of all the files in a particular folder.

First lets get a list of all files in our raster folder:

list.files(path="./data/raster/CWD/")
##  [1] "cwd2010apr.tfw"         "cwd2010apr.tif"         "cwd2010apr.tif.aux.xml"
##  [4] "cwd2010apr.tif.xml"     "cwd2010aug.tfw"         "cwd2010aug.tif"        
##  [7] "cwd2010aug.tif.aux.xml" "cwd2010aug.tif.xml"     "cwd2010dec.tfw"        
## [10] "cwd2010dec.tif"         "cwd2010dec.tif.aux.xml" "cwd2010dec.tif.xml"    
## [13] "cwd2010feb.tfw"         "cwd2010feb.tif"         "cwd2010feb.tif.aux.xml"
## [16] "cwd2010feb.tif.xml"     "cwd2010jan.tfw"         "cwd2010jan.tif"        
## [19] "cwd2010jan.tif.aux.xml" "cwd2010jan.tif.xml"     "cwd2010jul.tfw"        
## [22] "cwd2010jul.tif"         "cwd2010jul.tif.aux.xml" "cwd2010jul.tif.xml"    
## [25] "cwd2010jun.tfw"         "cwd2010jun.tif"         "cwd2010jun.tif.aux.xml"
## [28] "cwd2010jun.tif.xml"     "cwd2010mar.tfw"         "cwd2010mar.tif"        
## [31] "cwd2010may.tif"         "cwd2010nov.tif"         "cwd2010oct.tif"        
## [34] "cwd2010sep.tif"

Notice that .png and twf xml files are also listed. We can narrow this down by using a grep style pattern search.

Here we will only list files that have .tif in the name:

list.files(path="./data/raster/CWD/",
           pattern = '.tif')
##  [1] "cwd2010apr.tif"         "cwd2010apr.tif.aux.xml" "cwd2010apr.tif.xml"    
##  [4] "cwd2010aug.tif"         "cwd2010aug.tif.aux.xml" "cwd2010aug.tif.xml"    
##  [7] "cwd2010dec.tif"         "cwd2010dec.tif.aux.xml" "cwd2010dec.tif.xml"    
## [10] "cwd2010feb.tif"         "cwd2010feb.tif.aux.xml" "cwd2010feb.tif.xml"    
## [13] "cwd2010jan.tif"         "cwd2010jan.tif.aux.xml" "cwd2010jan.tif.xml"    
## [16] "cwd2010jul.tif"         "cwd2010jul.tif.aux.xml" "cwd2010jul.tif.xml"    
## [19] "cwd2010jun.tif"         "cwd2010jun.tif.aux.xml" "cwd2010jun.tif.xml"    
## [22] "cwd2010mar.tif"         "cwd2010may.tif"         "cwd2010nov.tif"        
## [25] "cwd2010oct.tif"         "cwd2010sep.tif"

almost there… we still have problems with .tif.* files like tif.xml or tif.aux.xml. We can limit these using what I will call an ANTI-wildcard. “$”. '.tif\$' will tell the computer the file needs to END in .tif. Let’s also store the path/file name by setting full.names=TRUE

list.files(path="./data/raster/CWD/",
           pattern = '.tif$', 
           full.names=TRUE)
##  [1] "./data/raster/CWD//cwd2010apr.tif" "./data/raster/CWD//cwd2010aug.tif"
##  [3] "./data/raster/CWD//cwd2010dec.tif" "./data/raster/CWD//cwd2010feb.tif"
##  [5] "./data/raster/CWD//cwd2010jan.tif" "./data/raster/CWD//cwd2010jul.tif"
##  [7] "./data/raster/CWD//cwd2010jun.tif" "./data/raster/CWD//cwd2010mar.tif"
##  [9] "./data/raster/CWD//cwd2010may.tif" "./data/raster/CWD//cwd2010nov.tif"
## [11] "./data/raster/CWD//cwd2010oct.tif" "./data/raster/CWD//cwd2010sep.tif"

Let’s store that information in a variable called CWDS:

CWDS = list.files(path="./data/raster/CWD/", 
                  pattern = 'tif$', 
                  full.names=TRUE)

Challenge: In teams read in and sum all of these rasters, then divide by the total number of files.

# Hint
jan = raster("./data/raster/CWD//cwd2010jan.tif")
feb = raster("./data/raster/CWD//cwd2010feb.tif")
# etc

# sum them then divide by total number of files  (need to include other months)
CWD_mn_2010 = (jan+feb ) / 2
### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

Finally, let’s write out our results:

# writeRaster(x=yourraster,filename="./output/CWD_mn_2010.tif")

Raster Stacks

R also has the capacity to handle multi-band or multi-layer raster ‘stacks’. raster stack In this case we will be using a stack to combine multiple files of the same type (CWD) over time (months of 2010). As you will see there are multiple advantages to handling data in stacks.

The only requirement is that all the images must have the same: - Image extent - Resolution - Projection

In this case our CWD data is already identical in shape and projection. All we need to do is use a list of files like CWDS that we created earlier.

CWDS
##  [1] "./data/raster/CWD//cwd2010apr.tif" "./data/raster/CWD//cwd2010aug.tif"
##  [3] "./data/raster/CWD//cwd2010dec.tif" "./data/raster/CWD//cwd2010feb.tif"
##  [5] "./data/raster/CWD//cwd2010jan.tif" "./data/raster/CWD//cwd2010jul.tif"
##  [7] "./data/raster/CWD//cwd2010jun.tif" "./data/raster/CWD//cwd2010mar.tif"
##  [9] "./data/raster/CWD//cwd2010may.tif" "./data/raster/CWD//cwd2010nov.tif"
## [11] "./data/raster/CWD//cwd2010oct.tif" "./data/raster/CWD//cwd2010sep.tif"

and create a stack from the data:

cwd_stack = stack(CWDS)

Now looking at the properties we can see that cwd_stack has 12 layers.

cwd_stack
## class      : RasterStack 
## dimensions : 1120, 872, 976640, 12  (nrow, ncol, ncell, nlayers)
## resolution : 1080, 1080  (x, y)
## extent     : -375034, 566726, -616284.9, 593315.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## names      : cwd2010apr, cwd2010aug, cwd2010dec, cwd2010feb, cwd2010jan, cwd2010jul, cwd2010jun, cwd2010mar, cwd2010may, cwd2010nov, cwd2010oct, cwd2010sep 
## min values :          0,          0,          0,          0,          0,          0,          0,          ?,          ?,          ?,          ?,          ? 
## max values :  150.96939,  205.35623,   47.32375,   65.93313,   32.79562,  233.96126,  219.87500,          ?,          ?,          ?,          ?,          ?

This allows us to do operations much easier. For instance we can just run arithmetic functions on the entire stack:

CWD_mn_2010 = mean(cwd_stack)
plot(CWD_mn_2010)

We can also do ‘global’ statistics like finding the ‘mean’ of all cells for any given month.

cellStats(cwd_stack, stat='mean', na.rm=TRUE)
##  cwd2010apr  cwd2010aug  cwd2010dec  cwd2010feb  cwd2010jan  cwd2010jul 
##  45.5497388 157.9214946   1.6478612   5.2903873   0.1809682 176.2727036 
##  cwd2010jun  cwd2010mar  cwd2010may  cwd2010nov  cwd2010oct  cwd2010sep 
## 143.6983408  27.6639933  90.9696350  19.2020629  27.7737959 120.8165948

Or other summary visualization like histograms on a month by month basis.

hist(cwd_stack)

Raster Vector Interactions

Raster cropping

Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors). Often the extent of input raster datasets is larger than the area of interest. In this case raster cropping and masking are useful for unifying the spatial extent of input data. Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.

We will use two objects to illustrate raster cropping:

  • A SpatRaster object srtm representing elevation (meters above sea level) in south-western Utah
  • A vector (sf) object zion representing Zion National Park

Both target and cropping objects must have the same projection. The following code chunk therefore not only reads the datasets from the spDataLarge package. This package needs to be installed directly from github using the following code:

# this has already been done for you, please don't run
#install.packages("spDataLarge", repos = "https://geocompr.r-universe.dev")

We will also need to install/import the sf package to handle vector data.

if(!require("spDataLarge")) install.packages("spDataLarge")
## Loading required package: spDataLarge
if(!require("sf")) install.packages("sf")

srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, crs(srtm))

Let’s take a look at the two files, note the use of add=TRUE in the second plot command. This adds the new object to the previous plot.

plot(srtm)
plot(zion, add=TRUE)

We use crop() from the raster package to crop the srtm raster. The function reduces the rectangular extent of the object passed to its first argument based on the extent of the object passed to its second argument.

srtm_cropped = crop(srtm, zion)
plot(srtm_cropped)
plot(zion, add=TRUE)

Note that the extent of the strm raster has been cropped to fit the zion polygon.

Related to crop() is the terra function mask(), which sets values outside of the bounds of the object are converted to NA, or a missing value. The following command therefore masks every cell outside of the Zion National Park boundaries:

srtm_masked = mask(srtm, zion)
plot(srtm_masked)

Importantly, we want to use both crop() and mask() together in most cases. This combination of functions would (a) limit the raster’s extent to our area of interest and then (b) replace all of the values outside of the area to NA.

srtm_cropped = crop(srtm, zion)
srtm_final = mask(srtm_cropped, zion)
plot(srtm_final)

Changing the settings of mask() yields different results. Setting updatevalue = 0, for example, will set all pixels outside the national park to 0. Setting inverse = TRUE will mask everything inside the bounds of the park (see ?mask for details)

srtm_inv_masked = mask(srtm, zion, inverse = TRUE)
plot(srtm_inv_masked)

Raster extraction

Raster extraction is the process of identifying and returning the values associated with a ‘target’ raster at specific locations, based on a (typically vector) geographic object location. The results depend on the type of selector used (points, lines or polygons) and arguments passed to the extract() function, which we use to demonstrate raster extraction. The reverse of raster extraction — assigning raster cell values based on vector objects — is rasterization.

The basic example is of extracting the value of a raster cell at specific points. For this purpose, we will use zion_points, which contain a sample of 30 locations within the Zion National Park.

# access some available data called zion_points
data("zion_points", package = "spDataLarge")
# now plot them all 
plot(srtm)
plot(zion, add=TRUE, col=NA)
plot(zion_points, add=TRUE, col='red')

The following command extracts elevation values from srtm and creates a data frame with points IDs (one value per vector’s row) and related srtm values for each point.

elevation = extract(srtm, zion_points)
elevation
##  [1] 1802 2433 1886 1370 1452 1635 1380 2032 1830 1860 1440 2145 1942 1691 1776
## [16] 2198 1820 1349 1758 1424 2159 1809 1826 1550 1799 2102 2118 1372 1905 1574

Now, we can add the resulting object to our zion_points dataset with the cbind() function:

zion_points = cbind(zion_points, elevation) # insert the elevation data into the orginal shapefile
zion_points
## Simple feature collection with 30 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2077 ymin: 37.16632 xmax: -112.8717 ymax: 37.43165
## Geodetic CRS:  WGS 84
## First 10 features:
##    elevation                   geometry
## 1       1802 POINT (-112.9159 37.20013)
## 2       2433 POINT (-113.0937 37.39263)
## 3       1886 POINT (-113.0246 37.33466)
## 4       1370 POINT (-112.9611 37.24326)
## 5       1452 POINT (-112.9898 37.20847)
## 6       1635 POINT (-112.8807 37.19319)
## 7       1380 POINT (-113.0505 37.24061)
## 8       2032 POINT (-113.0953 37.34965)
## 9       1830 POINT (-113.0362 37.31429)
## 10      1860 POINT (-113.2077 37.43165)

Note that now we have the elevation values for each one of the points associated with zion_points. For fun let’s use mapview to make an interactive map of both. Again we will use zcol to set the column used for symbology in the point map.

if(!require("mapview")) install.packages("mapview")
## Loading required package: mapview
mapview(srtm)+
mapview(zion_points, zcol='elevation')

Let’s tie this lesson back to our work with OSM. Here we are going to import the package osmdata to allow us to interact with the Overpass Turbo API. Let’s try to build a query using a bounding box from getbb(), where we search for Zion National Park, this bounding box then get passed to Overpass Turbo and restaurants will be added.

# import functions from osmdata library
if(!require("osmdata")) install.packages("osmdata")


# create a overpass turbo query
query = getbb(place_name="Zion National Park, Utah, USA") %>%
      opq() %>%
       add_osm_feature(key="amenity", value="restaurant")

# convert to a spatial data type sf 
restaurant = osmdata_sf(q=query)

Ok, bad request! Thats not too specific. Let’s try to see what happened. First let’s see if we were able to get a bounding box for Zion National Park.

getbb(place_name="Zion National Park, Utah, USA")

Looks like that is the problem. You can try to see if you can adjust the name to get a valid response. But its likely we will need to create a bounding box manually. Luckily we can get the bounding box from the srtm raster file using bbox(srtm), this then needs to be converted to a vector using c() in order to meet the requirement of the opq function. Look at ?opq, notice that we can pass a bounding box in the format bbox = c(xmin, ymin, xmax, ymax).

Note: bounding boxes can also be found using tools like bboxfinder.com.

Now let’s try again, using our new bounding box.

query =  opq(bbox = c(bbox(srtm))) %>%
       add_osm_feature(key="amenity", value="restaurant")

restaurant = osmdata_sf(q=query)
restaurant
## Object of class 'osmdata' with:
##                  $bbox : 37.1320834298579,-113.239583212784,37.5129167631658,-112.85208321281
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 32 points
##             $osm_lines : NULL
##          $osm_polygons : 'sf' Simple Features Collection with 4 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL

Let’s just focus on the points that we returned to us:

restaurant_points = restaurant$osm_points
restaurant_points
## Simple feature collection with 32 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2018 ymin: 37.16199 xmax: -112.8757 ymax: 37.27104
## Geodetic CRS:  WGS 84
## First 10 features:
##                osm_id                         name  addr.city addr.housenumber
## 880118222   880118222                         <NA>       <NA>             <NA>
## 883329842   883329842      Zion Pizza & Noodle Co.       <NA>             <NA>
## 904095151   904095151           Pioneer restaurant       <NA>             <NA>
## 937546158   937546158        Parallel Eighty-Eight       <NA>             <NA>
## 1579995974 1579995974                    Thai Sapa Springdale              198
## 2983983311 2983983311 Switchback Grille Steakhouse       <NA>             <NA>
## 2983983665 2983983665     Arkansas Al's Steakhouse       <NA>             <NA>
## 3641061615 3641061615                         <NA>       <NA>             <NA>
## 3641061616 3641061616                         <NA>       <NA>             <NA>
## 3641061617 3641061617                         <NA>       <NA>             <NA>
##            addr.postcode addr.state         addr.street addr.unit    amenity
## 880118222           <NA>       <NA>                <NA>      <NA>       <NA>
## 883329842           <NA>       <NA>                <NA>      <NA> restaurant
## 904095151           <NA>       <NA>                <NA>      <NA> restaurant
## 937546158           <NA>       <NA>                <NA>      <NA> restaurant
## 1579995974         84767         UT Zion Park Boulevard      <NA> restaurant
## 2983983311          <NA>       <NA>                <NA>      <NA> restaurant
## 2983983665          <NA>       <NA>                <NA>      <NA> restaurant
## 3641061615          <NA>       <NA>                <NA>      <NA>       <NA>
## 3641061616          <NA>       <NA>                <NA>      <NA>       <NA>
## 3641061617          <NA>       <NA>                <NA>      <NA>       <NA>
##                     cuisine level name.en name.fr
## 880118222              <NA>  <NA>    <NA>    <NA>
## 883329842       pasta;pizza     1    <NA>    <NA>
## 904095151              <NA>  <NA>    <NA>    <NA>
## 937546158              <NA>  <NA>    <NA>    <NA>
## 1579995974             thai  <NA>    <NA>    <NA>
## 2983983311 fish;steak_house  <NA>    <NA>    <NA>
## 2983983665         american  <NA>    <NA>    <NA>
## 3641061615             <NA>  <NA>    <NA>    <NA>
## 3641061616             <NA>  <NA>    <NA>    <NA>
## 3641061617             <NA>  <NA>    <NA>    <NA>
##                                                                                                    opening_hours
## 880118222                                                                                                   <NA>
## 883329842                                                                                      Mo-Su 16:00-21:00
## 904095151                                                                                                   <NA>
## 937546158                                                                                                   <NA>
## 1579995974                                                                                                  <NA>
## 2983983311 Apr-Oct 17:00-21:00; Dec-Feb Mo-We off; Dec-Feb Th-Su 17:00-20:00; Mar,Nov 17:00-20:00; Dec 24-25 off
## 2983983665                                                                                                  <NA>
## 3641061615                                                                                                  <NA>
## 3641061616                                                                                                  <NA>
## 3641061617                                                                                                  <NA>
##            outdoor_seating           phone takeaway
## 880118222             <NA>            <NA>     <NA>
## 883329842              yes +1-435-772-3815      yes
## 904095151             <NA>            <NA>     <NA>
## 937546158             <NA>            <NA>     <NA>
## 1579995974            <NA>    +14357720510     <NA>
## 2983983311            <NA> +1-435-772-3700     <NA>
## 2983983665            <NA>            <NA>     <NA>
## 3641061615            <NA>            <NA>     <NA>
## 3641061616            <NA>            <NA>     <NA>
## 3641061617            <NA>            <NA>     <NA>
##                                        website                   geometry
## 880118222                                 <NA> POINT (-113.0003 37.18914)
## 883329842      http://www.zionpizzanoodle.com/ POINT (-112.9984 37.18924)
## 904095151                                 <NA> POINT (-112.9976 37.18926)
## 937546158                                 <NA> POINT (-113.0062 37.17994)
## 1579995974        http://www.thaisapazion.com/ POINT (-112.9918 37.19809)
## 2983983311 http://switchbackgrille.com/grille/ POINT (-113.0026 37.18402)
## 2983983665                                <NA>  POINT (-113.014 37.16757)
## 3641061615                                <NA> POINT (-112.8759 37.27104)
## 3641061616                                <NA>  POINT (-112.8757 37.2709)
## 3641061617                                <NA> POINT (-112.8758 37.27078)

That is too much data. Let’s select only those columns we are interested in:

restaurant_points = subset(restaurant_points, select = c(name,phone ))
restaurant_points
## Simple feature collection with 32 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2018 ymin: 37.16199 xmax: -112.8757 ymax: 37.27104
## Geodetic CRS:  WGS 84
## First 10 features:
##                                    name           phone
## 880118222                          <NA>            <NA>
## 883329842       Zion Pizza & Noodle Co. +1-435-772-3815
## 904095151            Pioneer restaurant            <NA>
## 937546158         Parallel Eighty-Eight            <NA>
## 1579995974                    Thai Sapa    +14357720510
## 2983983311 Switchback Grille Steakhouse +1-435-772-3700
## 2983983665     Arkansas Al's Steakhouse            <NA>
## 3641061615                         <NA>            <NA>
## 3641061616                         <NA>            <NA>
## 3641061617                         <NA>            <NA>
##                              geometry
## 880118222  POINT (-113.0003 37.18914)
## 883329842  POINT (-112.9984 37.18924)
## 904095151  POINT (-112.9976 37.18926)
## 937546158  POINT (-113.0062 37.17994)
## 1579995974 POINT (-112.9918 37.19809)
## 2983983311 POINT (-113.0026 37.18402)
## 2983983665  POINT (-113.014 37.16757)
## 3641061615 POINT (-112.8759 37.27104)
## 3641061616  POINT (-112.8757 37.2709)
## 3641061617 POINT (-112.8758 37.27078)

Now let’s plot it out and assign the restaurant names as labels:

mapview(srtm)+
mapview(restaurant_points,
        label = restaurant_points$name,
        legend=F)

Let’s extract the elevation of each restaurant, following the same code as earlier:

elevation = extract(srtm, restaurant_points)
elevation
##  [1] 1197 1194 1194 1173 1212 1184 1175 1984 1984 1984 1984 1153 1153 1162 1162
## [16] 1210 1202 1200 1200 1202 1195 1194 1197 1184 1197 1197 1200 1089 1318 1192
## [31] 1191 1196

Now, we can add the resulting list of elevation values to our restaurant_points dataset with the cbind() function:

restaurant_points = cbind(restaurant_points, elevation) # insert the elevation data into the orginal shapefile
restaurant_points
## Simple feature collection with 32 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2018 ymin: 37.16199 xmax: -112.8757 ymax: 37.27104
## Geodetic CRS:  WGS 84
## First 10 features:
##                                    name           phone elevation
## 880118222                          <NA>            <NA>      1197
## 883329842       Zion Pizza & Noodle Co. +1-435-772-3815      1194
## 904095151            Pioneer restaurant            <NA>      1194
## 937546158         Parallel Eighty-Eight            <NA>      1173
## 1579995974                    Thai Sapa    +14357720510      1212
## 2983983311 Switchback Grille Steakhouse +1-435-772-3700      1184
## 2983983665     Arkansas Al's Steakhouse            <NA>      1175
## 3641061615                         <NA>            <NA>      1984
## 3641061616                         <NA>            <NA>      1984
## 3641061617                         <NA>            <NA>      1984
##                              geometry
## 880118222  POINT (-113.0003 37.18914)
## 883329842  POINT (-112.9984 37.18924)
## 904095151  POINT (-112.9976 37.18926)
## 937546158  POINT (-113.0062 37.17994)
## 1579995974 POINT (-112.9918 37.19809)
## 2983983311 POINT (-113.0026 37.18402)
## 2983983665  POINT (-113.014 37.16757)
## 3641061615 POINT (-112.8759 37.27104)
## 3641061616  POINT (-112.8757 37.2709)
## 3641061617 POINT (-112.8758 37.27078)

Extract Values Along a Line

Raster extraction also works with line selectors. Then, it extracts one value for each raster cell touched by a line. In this case, the best approach is to split the line into many points and then extract the values for these points. To demonstrate this, the code below creates zion_transect, a straight line going from northwest to southeast of the Zion National Park

Let’s create a line based on two points c(-113.2, -112.9) and c(37.45, 37.2) which are xy pairs in lat lon. Note that we are assigning a projection using an EPSG code of 4326. EPSG codes are unique identifying codes that are available for all common projections. To look up EPSG codes please go to https://epsg.io.

zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
  st_linestring() %>%
  st_sfc(crs =  'EPSG:4326')%>%
  st_sf()
zion_transect
## Simple feature collection with 1 feature and 0 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS:  WGS 84
##                         geometry
## 1 LINESTRING (-113.2 37.45, -...

Now let’s see if its showing up in the right place

mapview(srtm)+
mapview(zion_transect, 
        color='white')

The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike. The method demonstrated below provides an ‘elevation profile’ of the route (the line does not need to be straight), useful for estimating how long it will take due to long climbs.

Because our line is made up of only two points and extract works best with points, we need to break our longer line into lots of points.

The first step is to add a unique id for each transect. Next, with the st_segmentize() function we can add points along our line(s) with a provided density (dfMaxLength) and convert them into points with st_cast().

Let’s look at st_segmentize help to see what units dfMaxLength is in:

?st_segmentize

Now lets create our collection of points:

zion_transect$id = 1:nrow(zion_transect)
zion_transect = st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect = st_cast(zion_transect, "POINT")
zion_transect
## Simple feature collection with 257 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS:  WGS 84
## First 10 features:
##     id                   geometry
## 1    1       POINT (-113.2 37.45)
## 1.1  1 POINT (-113.1988 37.44902)
## 1.2  1 POINT (-113.1976 37.44805)
## 1.3  1 POINT (-113.1965 37.44707)
## 1.4  1  POINT (-113.1953 37.4461)
## 1.5  1 POINT (-113.1941 37.44512)
## 1.6  1 POINT (-113.1929 37.44415)
## 1.7  1 POINT (-113.1918 37.44317)
## 1.8  1  POINT (-113.1906 37.4422)
## 1.9  1 POINT (-113.1894 37.44122)

Now, we have a large set of points, and we want to derive a distance between the first point in our transects and each of the subsequent points. In this case, we only have one transect, but the code, in principle, should work on any number of transects:

zion_transect = zion_transect %>%
  group_by(id) %>%
  mutate(dist = st_distance(geometry)[, 1]) 
zion_transect
## Simple feature collection with 257 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS:  WGS 84
## # A tibble: 257 × 3
## # Groups:   id [1]
##       id             geometry  dist
##  * <int>          <POINT [°]>   [m]
##  1     1       (-113.2 37.45)    0 
##  2     1 (-113.1988 37.44902)  150.
##  3     1 (-113.1976 37.44805)  300.
##  4     1 (-113.1965 37.44707)  450.
##  5     1  (-113.1953 37.4461)  600.
##  6     1 (-113.1941 37.44512)  750.
##  7     1 (-113.1929 37.44415)  901.
##  8     1 (-113.1918 37.44317) 1051.
##  9     1  (-113.1906 37.4422) 1201.
## 10     1 (-113.1894 37.44122) 1351.
## # … with 247 more rows

Finally, we can extract elevation values for each point in our transects and combine this information with our main object.

zion_elev = extract(srtm, zion_transect)
zion_transect = cbind(zion_transect, zion_elev)
zion_transect
## Simple feature collection with 257 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS:  WGS 84
## First 10 features:
##    id          dist zion_elev                   geometry
## 1   1    0.0000 [m]      2001       POINT (-113.2 37.45)
## 2   1  150.0961 [m]      2033 POINT (-113.1988 37.44902)
## 3   1  300.1922 [m]      2020 POINT (-113.1976 37.44805)
## 4   1  450.2883 [m]      1989 POINT (-113.1965 37.44707)
## 5   1  600.3844 [m]      1902  POINT (-113.1953 37.4461)
## 6   1  750.4805 [m]      1869 POINT (-113.1941 37.44512)
## 7   1  900.5766 [m]      1827 POINT (-113.1929 37.44415)
## 8   1 1050.6728 [m]      1767 POINT (-113.1918 37.44317)
## 9   1 1200.7689 [m]      1772  POINT (-113.1906 37.4422)
## 10  1 1350.8650 [m]      1722 POINT (-113.1894 37.44122)

We can now see that we have both the distance from the first point in our line, and the elevation of each point along that line.

Let’s plot the transect and see what it looks like:

plot(x=zion_transect$dist, 
     y = zion_transect$zion_elev, 
     type='l', 
     col='darkred',
     main='Elevation Profile',
     xlab = 'Distance',
     ylab = 'Elevation (m)')

Step 1 - Create a query to pull all the ‘primary’ roads in Zion National Park and save them in a variable called roads. Make sure to look at the tagging convention on OSM’s map features directory. Also you need to set the bbox from the query, you will want to do it based on the srtm raster as we did earlier.

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# # hint
# query =  opq(bbox =  ?????) %>%
#        add_osm_feature(key="?????", value="?????")
# 
# roads = osmdata_sf(q=query)

Step 2: osmdata_sf returns points (roads$osm_points), lines (roads$osm_lines), and polygons(roads$osm_polygons). Create a new variable called primary_roads that holds the road line features.

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# print(primary_roads)

Step 3: Use the filter() function from dpylr to isolate only roads with the name Zion Park Boulevard, and store the output in a variable called park_road:

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# print(park_road)

Let’s take a look at the data:

park_road

Notice that there are multiple segments making up this one road. Roads are typically broken into segments to allow for connectivity between other roads, and to maintain data like speed limits - which likely vary along the road.

Step 4: We need to ‘dissolve’ all of these road features into a single line representing ‘Zion Park Boulevard’. In order to do this we are going to tell the computer to create a new geometry column using summarize, based on the st_union of the multiple geometry rows. This is one that we haven’t talked about yet, so I will just do it for you.

park_road = park_road %>% 
            summarize(geometry = st_union(geometry))
plot(park_road)

print(park_road)
## Simple feature collection with 1 feature and 0 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -113.0368 ymin: 37.16021 xmax: -112.9902 ymax: 37.19982
## Geodetic CRS:  WGS 84
##                         geometry
## 1 LINESTRING (-113.0368 37.16...

Notice there is now only one row of data and our multiple road segments are now the union.

Step 5: Now we can create our column called id, then use st_segmentize to break the road into multiple points at a fixed distance (250 or 500 works well), then we need to cast that point object back into an sf point object using st_cast.

See if you can find the correct chunk of code

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# park_road$id = 1:nrow(???)
# park_road = st_segmentize(park_road, dfMaxLength = ????)
# park_road = st_cast(park_road, "?????")
# park_road

Step 6: Let’s calculate the distances from the first point to all following points:

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# park_road = park_road %>%
#   group_by(?????) %>%
#   mutate(dist = st_distance(?????)[, 1]) 
# park_road

Step 7: Now we can use the extract function to grab the values from the raster srtm at each one of the points making up park_road, then we need to cbind those values back into the park_road_elev dataset.

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# park_road_elev = extract(????, park_road)
# park_road_elev = cbind(????, ????)
# park_road_elev

Step 8: Plot the park_road_elev$dist on the x axis and park_road_elev$park_road_elev on the y axis, and set the type to 'line'.

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# plot(x = park_road_elev$?????, 
#      y = park_road_elev$?????, 
#      type='line')